Improved gross primary production estimation in rice fields through integrated multi‐scale methodologies

Abstract Understanding productivity in agricultural ecosystems is important, as it plays a significant role in modifying regional carbon balances and capturing carbon in the form of agricultural yield. This study in particular combines information from flux determinations using the eddy covariance (EC) methodology, process‐based modeling of carbon gain, remotely (satellite) sensed vegetation indices (VIs), and field surveys to assess the gross primary production (GPP) of rice, which is a primary food crop worldwide. This study relates two major variables determining GPP. The first is leaf area index (LAI) and carboxylation capacity of the rice canopy (Vcuptake), and the second being MODIS remotely sensed vegetation indices (VIs). Success in applying such derived relationships has allowed GPP to be remotely determined over the seasonal course of rice development. The relationship to VIs of both LAI and Vcuptake was analyzed first by using the regression approaches commonly applied in remote sensing studies. However, the resultant GPP estimations derived from these generic models were not consistently accurate and led to a large proportion of underestimations. The new, alternative approach developed to estimate LAI and Vcuptake uses consistent development curves for rice (i.e., relies on consistent biological regulations of plant development). The modeled GPP based on this consistent development curve for both LAI and Vcuptake agreed with R 2 from 0.76 to 0.92 (within the 95% confidence interval). The results of this study demonstrate that improved linkages between ground‐based survey data, eddy flux measurements, process‐based models, and remote sensing data can be constructed to estimate GPP in rice paddies. This study suggests further that the conceptual application of the consistent development curve, such as the combining of different scale measurements, has the potential to predict GPP better than the common practice of utilizing simple linear models, when seeking to estimate the critical parameters that influence carbon gain and agricultural yields.


| INTRODUC TI ON
Gross primary production (GPP) represents the total carbon assimilation by vegetation through photosynthesis, playing a vital role in crop yield, carbon balance, and ecosystem functioning (Beer et al., 2010;Moureaux et al., 2008;O'Sullivan et al., 2020). GPP refers to the total uptake of carbon dioxide in photosynthesis or CO 2 assimilation, which drives subsequent ecosystem processes, including plant growth and agricultural yield (Schulze, 2006). Accurate estimation of GPP is essential for understanding the impacts of climate variability on agricultural productivity and ensuring global food security (Huntzinger et al., 2017;Liu et al., 2020). Rice (Oryza sativa L.) is a staple food crop for over 50% of the world population and occupies approximately 150 million hectares of land, making the accurate estimation of rice GPP critical Ray et al., 2013).
In recent years, significant advancements have been made in estimating GPP using eddy covariance (EC) methods, process-based models, and satellite remote sensing (Baldocchi, 2014;Li et al., 2018;Peng & Gitesol, 2011Zhang et al., 2018). The increased availability of high-resolution satellite data, such as Landsat-8 and Sentinel-2, as well as rapid developments in UAV technology, has improved the spatial resolution of GPP estimations (Chaves et al., 2020;Wulder et al., 2019). However, estimating GPP in agroecosystems like rice paddies remains a challenge due to the incompatibility of spatial scales of flux measurements and satellite imagery resolution, as well as the influence of local climate, water availability, and varied management practices (Wattenbach et al., 2010;Yu et al., 2017).
Recent studies have reported the integration of advanced remote sensing techniques, such as near-infrared vegetation indices (NIRv) and near-infrared vegetation index with photosynthesis (NIRvP), to better capture rapid changes in crop development and improve GPP estimation (Badgley et al., 2017;Luo et al., 2019). The main goal of this study is to assess GPP at rice paddies in Asia and Europe, integrating recent advancements in remote sensing, process-based modeling, and eddy covariance measurements to provide a more accurate estimation of rice GPP (Bouman et al., 2007;Gan et al., 2021).
While recent research has made significant advancements in GPP monitoring and modeling, the current study differentiates itself by integrating multiple data sources and techniques for estimating GPP in rice paddies, addressing the unique challenges posed by this agroecosystem. By combining ground-based survey data, eddy flux measurements, process-based models, and remote sensing data, this study offers a more accurate and comprehensive approach for GPP estimation in rice paddies (Ogutu et al., 2013;Peng & Gitesol, 2012).
This integrative approach not only enhances the understanding of the impacts of climate variability on agricultural productivity but also provides valuable insights for informed crop management decisions and ensuring global food security (Revill et al., 2013).
In this study, we use the physiologically based PIXGRO canopy model, combined with high-resolution remote sensing data and eddy covariance measurements, to estimate rice paddy GPP (Li et al., 2018;Tenhunen et al., 2009). We aim to establish new relationships between remotely sensed vegetation indices, such as NDVI, and GPP and biomass, enabling improved parameterization of processbased models at the landscape scale (Gamon et al., 1995;Nemani et al., 2003;Wang et al., 2004).
The objectives of this study are: 1. To use NDVI to examine rice paddy phenology and establish relationships between observed LAI and vegetation indices (Lee et al., 2016) 2. To estimate the critical PIXGRO parameter Vc uptake of rice paddies from EC daily values in conjunction with leaf area index (LAI), which controls for canopy carbon fixation in response to multiple environmental factors (Owen et al., 2007) 3. To define the relationship between the daily NDVI and modeling parameter Vc uptake , establishing best-fit models to estimate GPP in the absence of observations of CO 2 exchange and LAI.
The process-based model focused on the physiological processbased canopy sub-model of PIXGRO (Tenhunen et al., 2009) that links flux observations from the eddy covariance studies with ecosystem physiology represented as the capacity for CO 2 exchange (i.e., canopy content and activity of Rubisco) and plant phenology (i.e., leaf area index). PIXGRO is designed as a tool for bridging between measured gas exchange fluxes, derived parameters for carboxylation capacity, seasonal changes in biomass and structure in the case of herbaceous and crop plants, and crop yields, considering specific ecophysiological behavior of individual species (Adiku et al., 2006). The canopy sub-model of PIXGRO (hereafter referred to as "canopy model") calculates the dynamics of whole-ecosystem CO 2 and H 2 O exchange (Reichstein, 2001).
We hypothesize that combining the physiologically based PIXGRO canopy model with GPP obtained from EC rice paddy sites allows for a more accurate definition of the seasonal course of the critical parameter Vc uptake . The seasonal course of the parameter Vc uptake , together with meteorological data, enables a more efficient reproduction of observed data that then helps more clearly define GPP. We also hypothesize that Vc uptake , estimated using seasonally observed LAI, performs better in estimating GPP than Vc uptake using a constant LAI. Estimating GPP through this best-fit model for Vc uptake and LAI in dependence on NDVI subsequently provides more accurate GPP predictions (Li et al., 2008(Li et al., , 2018Zhang et al., 2018).

| Study sites
The analysis was conducted using available rice paddy data from South Korea (Haean Catchment) and Japan (Mase) in Asia and Spain  Table 1. The detailed information on landscape structure is described in Seo et al. (2014).
The total area of the catchment is 64 km 2 , consisting of 58% forested mountain area, 30% agricultural area, and 12% as residential, riparian, field margins, and farm road area according to land surveys (Arnhold et al., 2014). The agricultural area is characterized as a mosaic patchwork of fields, with a dominance of dry-land fields (22% of the total area) and rice paddy fields (8%) as the remaining area. Rice paddies (Oryza sativa L., cv. Odae) are cultivated at <500 m a.s.l. in the catchment (Choi et al., 2010).

| Mase, Japan (MSE)
The Mase site is located in the rural area (36° 03′ 14′′ N, 140° 01′ 38″ E, 15 m a.s.l.) of Tsukuba City in Central Japan. The rice paddy (Oryza sativa L.) was ca. 2 km 2 and was managed as a single rice-cropping field following practices common in the area (Saito et al., 2005). In this study, the EC flux and meteorological data from 2002 to 2005 obtained from AsiaFlux (https://db.cger.nies.go.jp/asiaf luxdb/) were included.

| El Saler-Sueca, Spain (ESES2)
The El Saler Sueca site is located in the protected wetland area of La Albufera Natural Park in the Valencia region of Spain (39° 16′ 32′′ N, 0° 18′ 55″ E, 10 m a.s.l.). El Saler, which is in a sub-arid Mediterranean climate, experiences hot summers with almost no rain and cold winters with substantial rainfall. The rice paddy was ca. 15 km 2 and had been managed in same way for 200 years Moors et al., 2010 Note: Total Rg is total global radiation (in parentheses, standard deviation), mean Ta is mean air temperature (in parentheses, standard deviation), total P is total rainfall, planting, and harvest date (day of year-DOY), maximum LAI, and DOY of maximum LAI. Meteorological condition considers the period of crop growth (from the planting to the harvest). (Reichstein et al., 2003;Tenhunen et al., 1994;Wang et al., 2003).
The model is single-layered model. It estimates light interception and CO 2 exchange rates of canopy foliage for sun and shaded light classes half-hourly, which is then compared to EC measurements (Owen et al., 2007). The model is driven by meteorological data, for example, global radiation (R g ), air temperature (T a ), vapor pressure deficit (VPD), wind speed, air pressure, and atmospheric CO 2 concentration, and requires estimated values for LAI. Total shortwave radiation on the sunlit leaves is the sum of direct, sky diffuse, and multiple scattered radiation, whereas on the shaded leaves, it is only the sum of sky diffuse and multiple scattered radiation (see Equations 2-5 in Owen et al., 2007). The foliage orientation function (G) was set at 0.5 and the influence of clumping (Ω) at 0.9 for croplands in these equations. In order to account for the effect of the canopy on light interception, we expanded LAI to plant area index (PAI), which is the sum of LAI and stem area index (SAI) (i.e., PAI = LAI + SAI). SAI of the crop is calculated as 14% of LAI, whereas SAI of the rice is set at 0.01 (see details in Owen et al., 2007).

Simulation of gross photosynthesis follows Farquhar and
Caemmerer (1982), as modified for practical field applications by Harley and Sharkey (1991). It is based on Ribulose-1,5bisphosphate-carboxylase-oxygenase (Rubisco) enzyme reactions, where the rate of CO 2 fixation is limited by either the regeneration of Ribulose-1,5-biphosphate (RuBP) at low light intensity and/or high internal CO 2 concentration or by Rubisco activity and CO 2 /O 2 concentration at saturated light and low internal CO 2 concentration. The key parameter of the model is Rubisco maximum carboxylation rate (Vc max ) at 25°C, while all temperature dependencies are fixed in relation to this rate. When comparing predicted GPP to EC measured values, a best fit is obtained for this key parameter. RuBP reduction capacity, dark respiration capacity, and light utilization efficiency of the canopy are assumed to be proportional to Vc max . Given that fixed temperature dependencies and process proportionalities are used and that assumptions are made about canopy structure and light interceptions, a lumped parameter, Vc uptake , that is assumed to control overall carbon fixation rather than direct enzyme related parameter Vc max is obtained from the statistical fitting procedure.
The model formulation follows Farquhar and Caemmerer (1982) further; net photosynthesis (P net ) is obtained using where Γ * is CO 2 compensation point in the absence of mitochondrial respiration, w c is the carboxylation rate supported by Rubisco enzyme, w j is the carboxylation rate supported by the actual electron transport rate, and R d is the respiration occurring in mitochondria without light.
c i is the internal CO 2 concentration based on Fick's Law for molecular diffusion of CO 2 through the stomata and boundary layer and is calculated from the following equation.
where c s is the CO 2 concentration at the surface of the leaf and g s is the stomatal conductance according to modified Ball-Berry equation (Ball et al., 1987;Harley & Sharkey, 1991).
where g s,min is the minimum stomatal conductance, rH is relative humidity, and gfac is a constant representing stomatal sensitivity in relation to CO 2 assimilation. It has been evaluated for different species from chamber experiments (Sala & Tenhunen, 1994Tenhunen et al., 1990).
w c is the carboxylation rate supported by Rubisco enzyme, calculated as: where Vc is the maximum rate of carboxylation, K c is the Michaelis con- where Vc max is the maximum rate of carboxylation capacity at 25°C, H a is the activation enthalpy of carboxylation, T k is the estimated the leaf temperature in the current model iteration step, R is the universal gas constant, S is an enthalpy term for deactivation, and H a is the deactivation enthalpy of carboxylation. w j is calculated as: where P m is the maximum potential rate of RuBP production. P m is calculated following the Smith equation (cf. Tenhunen et al., 1976): where alpha is the average leaf light utilization efficiency without photorespiration, I is the incident PPFD, and P ml is the CO 2 and light saturated temperature-dependent potential RuBP regeneration rate as described in Falge et al. (1997).
As indicated above, rather than using the notation of Vc max at the leaf level, Vc uptake is referred to as the estimate of the maximum rate of carboxylation capacity at the canopy level when fitting the model to EC data. Estimated Vc uptake is extracted by minimizing the sum of residual least squares, that is, the Levenberg-Marquardt algorithms of the PV-WAVE routine in the comparison of model predictions with EC observations at half hour intervals. The light utilization efficiency, alpha, can also be estimated from the EC data as an additional fitting parameter. In this study, however, we assumed alpha to be proportional to Vc uptake following Owen et al. (2007). The relationship of alpha and Vc uptake is Vc uptake estimation was carried out on a daily basis over the course of the growing season for each crop data set. We determined that the predicted Vc uptake was occasionally unrealistically high during the early growing season (>200 μmol m-2 leaf area s-1) when LAI was low (<1), but GPP was positive (>50 μmol m-2 leaf area s-1).
This results from errors in the flux determination or in LAI estimation, values that are critical in the analysis of CO 2 uptake (it may occur, for example, if the flux footprint is different from the area used in LAI estimation). In order to eliminate artificially high estimates of Vc uptake , the current work includes only flux data where LAI is >1 in the model fitting procedure. More information on the physiologically based model can be found in Owen et al. (2007).

| Physiological parameter estimates
The leaf physiological parameters applied as constants and those controlling temperature dependencies were obtained in previous studies on leaf physiology. These values are identified and shown in Table 1 (from Falge et al., 1996;Harley & Sharkey, 1991;Sala & Tenhunen, 1996;Tenhunen et al., 1990). These parameters describe temperature and light dependencies and response of stomata (Tenhunen et al., 1990). Gross primary production (GPP) is calculated with the single-layered canopy sub-model as described in section 2.2.1. Inputs to the model are half-hourly global radiation, air temperature, relative humidity, and precipitation. Matrices for all meteorological drivers are prepared previous to analysis runs (estimated in separate routines and stored outside of the model) and are input to the model according to the half-hourly simulation time step. Of primary concern in this study was estimation of a single key parameter, Vc uptake , that describes change in overall canopy CO 2 uptake capacity.
Seasonal variation of Vc uptake was estimated via model fitting with the detailed canopy sub-model and EC GPP data (Owen et al., 2007).
The method was performed using the functions NLINLSA and Consistent seasonal trends in the key physiological parameter describing CO 2 uptake capacity, Vc uptake , are found for functional crop types, for example, root crops and rice as a grain crop , which aids in parameterization according to the land use.

| Leaf area index estimates
Leaf area index (LAI) is the key parameter related to plant phenology and measured during the growing season at the rice paddy sites, Haean, Mase, and El Saler. Since measured LAI has in previous studies been shown to be correlated with spectral reflectance (Fan et al., 2008;Pontailler et al., 2003;Stenberg et al., 2004;Xiao et al., 2002), NDVI was applied to estimate LAI. The relationship between vegetation indices and LAI is generally applicable or

| Physiological parameter estimation
The  low RMSE (below 9.93), and low CV (3%-25%) reflect the agreement illustrated in Figure 2 and further demonstrate that LAI obtained by the consistent development method is useful in the estimation of GPP at the study sites.

| GPP estimation with the best-fit model
Daily GPP as reproduced by the model, which was LAI and Vc uptake both dependent on consistent development curves, is in Figure 3.
Accumulated observed GPP during the growing season varied from 672 gC m −2 d −1 to 1294 gC m −2 d −1 , although the period for comparison varied slightly. Carbon uptake appeared to increase from Korea, to Japan, and to Spain due to different climate conditions.
Modeled values ranged from 670 gC m −2 d −1 to 1020 gC m −2 d −1 with R 2 above 0.79, RMSE below 3.48, CV above 38.02% and modeling efficiency (MF) between 0.92 and 0.72 (Table 3). Simulated GPP was in general under-estimated as expected due to the remaining difficulties in estimating Vc uptake in dependence on NDVI.
Deviations from observation over the course of the season at each site reflect the differences found for observed and predicted Vc uptake as represented in Figure 3 or in Table 3. Although the alternative method based on consistent phenological development did not lead to an improvement in GPP prediction, further development of the response curve shown in Figure 2 may lead to better success as discussed later.

| DISCUSS ION
In developing a methodology for estimating GPP at a large scale, it was determined that inputting remotely sensed information into vegetation canopy models would help to achieve more accurate estimates. However, the reflectance observed by remote sensing is influenced by both physiological parameters and vegetation phenological development (i.e., changing canopy structure). The relative importance of these in determining vegetation indices is largely unknown. While in this study separation of their effects was not attempted, an approach of consistent development curves is used to relate NDVI stepwise both to LAI change and to physiological change over the course of the season. Seasonal changes in LAI and in the key physiological model parameter Vc uptake , which reflects canopy carboxylation capacity at a given LAI (described as dependent on the vegetation index NDVI). The relative successes and failures in these efforts are discussed below.

| Use of NDVI to estimate Vc uptake for rice paddies
The eddy covariance methodology (EC) applied at rice paddy sites quantifies seasonal changes in GPP in relation to prevailing meteorological conditions (Kwon et al., 2010). By fitting the PIXGRO canopy model routine with flux data and with known LAI and meteorological conditions on a daily basis, an estimated seasonal time course for the key physiological parameter Vc uptake (i.e., canopy carboxylation capacity) is obtained. Owen et al. (2007) obtained such estimates for many EC monitoring sites by assuming a constant LAI at the observed maximum level. They carried out their analysis with constant LAI because measured time courses TA B L E 2 Statistics for the linear correlation between Vc uptake_org and Vc uptake . a is a slope, b (μmol CO 2 m −2 s −1 ) is an intercept, R 2 is the coefficient of determination, RMSE (μmol CO 2 m −2 s −1 ) is root mean square error, and CV is coefficient of variation. F I G U R E 2 A scaled general seasonal curve for Vc uptake utilizing data from all rice paddy sites. The scaling results by dividing estimated Vc uptake by maximum Vc uptake for each site and year. The black closed circle indicates the relationship between the scaled Vc uptake and DOY based on time shifts according to maximum NDVI. The solid line indicates the general seasonal curve determined by generalized addictive model (GAM).
for LAI (periodic measurements) did not exist at many of the study sites. While an assumption of a constant LAI at maximum level may allow a more complete use of data from EC study sites and in a way lead to reasonable spatial models for GPP, the objective of the current study is to use remote sensing first to determine structural change (the seasonal course of LAI). This requires comparing

| CON CLUS ION
This study provides insights with respect to achieving improved methods for assessing GPP in rice paddy ecosystems through the use of process-based models and remote sensing. The MODIS product provided by NASA underestimates and poorly represents GPP, especially in croplands Yan et al., 2009;Zhang et al., 2008 (Table 3), although limited so far to periods where LAI >1.0.
Photosynthetic capacity in relation to VIs changes over the course of a growing season together with the carbon-to-nitrogen ratio of leaves in the plant canopy. Nitrogen in foliage determines vegetation's investments into chlorophyll and rubisco (Bonan, 2015), such as machinery for photosynthetic processes. This may decrease over the course of the season with little influence on reflectance properties. Where data were available, this study found that the carbonto-nitrogen ratios, often accompanied changes in vegetation indices.
These types of physiological changes (reflecting the relationship between Vc uptake and NDVI) may account for the hysteresis effects that were observed in particular data sets.
Improvements in the modeling approach require that more detailed biological supplementary information should accompany EC studies in the future (e.g., frequent spatial samplings of crop canopy structures, local on-the-ground measurements of VIs, monitoring of physiological characteristics such as C/N ratios, chlorophyll content). On the one hand, detailed studies must be designed to support an extension of the approach described here, moving away from empiricism and linking remote sensing to the biological processes being observed. Ultimately, however, the estimation of crop GPP should be possible and should be undertaken without the need for additional field surveys, since these cannot be conducted on a large scale.

ACK N OWLED G M ENTS
We would like to express our heartfelt gratitude to the Warm- Complex Terrain and Ecological Heterogeneity (GRK 1565/1), at the University of Bayreuth, Germany. Our gratitude extends to our colleagues and fellow researchers for their insightful inputs and discussions that greatly enriched this work. Finally, our thanks go to the reviewers and editors for their constructive comments and suggestions, which have substantially improved the quality of our manuscript.

FU N D I N G I N FO R M ATI O N
Funding was provided by National Institute of Forest Science (FE100-2022-04-2023).

CO N FLI C T O F I NTE R E S T S TATE M E NT
The authors declare no conflict of interest.

DATA AVA I L A B I L I T Y S TAT E M E N T
Data derived from public domain resources.